Contribution of the Large Magellanic Cloud 
to the Galactic Warp 



Toshio Tsuchiya 

Astronomisches Rechen- Institute, Monchhofstrafie 12-14-, D69120 Heidelberg, 

Germany 



Abstract 

Multi-scale interaction between the LMC, the Galactic halo, and the disk is exam- 
ined with N-body simulations, and precise amplitudes of the Galactic warp excita- 
tion are obtained. The Galactic models are constructed most realistically to satisfy 
available observational constraints on the local circular velocity, the mass, surface 
density and thickness of the disk, the mass and size of the bulge, the local density of 
the halo matter at the solar radius, and the mass and orbit of the LMC. The mass 
of the halo within i? = 50 kpc is set to about 5 x lO^^M©. Since the observational 
estimate of the mass distributed in outer region has large ambiguity, two extreme 
cases are examined; M(< 170kpc) = 2.1 and 0.9 x lO^^M©. LMC is orbiting in a 
ellipse with apocentric radii of 100 kpc, thus the main difference between our two 
models is the mass density in the satellite orbiting region, so that our study can 
clarify the role of the halo on excitation of the warp. 

By using hybrid algorithm (SCF-TREE) I have succeeded to follow the evolution 
with millions of particles. The orbiting satellite excites density enhancement as a 
wake, and the wake exerts a tidal force on the disk. Because of the additional torque 
from the wakes in the halo, the amplitudes of the induced warps are much larger 
than the classical estimate by Hunter & Toomre (1969), who considered only the 
direct torque from the LMC. The obtained amplitudes of m = 0, 1, 2 warps in the 
larger halo model show very good agreement with the observed amplitude in the 
Milky Way. This result revives the LMC as a possible candidate of the origin of the 
Galactic warp. Our smaller halo model, however, yield only weak warps in all the 
harmonic modes. Therefore the halo still has significant influence on excitation of 
warp even in the interaction scenario for excitation of warps. 
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1 Introduction 



The origin of galactic disk warping, which is typically vertical bending of disks 
in a integral shape, is a long-standing problem. Even after 40 year efforts 
to explain the phenomena since the discovery of the warp in the Milky Way 
(Burke 1957, Kerr 1957), there is still no generally acceptable theory. A general 
review of the problem can be found in Binney (1992). 

There are several observational facts that any convincing theory must explain. 
The first point is the longevity of the warps. Observations of neutral hydrogen 
layers (Bosma 1991) and stellar disks (Sanchez-Saavedra et al. 1990, Reshet- 
nikov & Combes 1998) have revealed that about a half of spiral galaxies ex- 
hibit warps. Such high frequency of warps indicate that the warps are long- 
lived with about the same age as galaxies, otherwise some mechanism keeps 
exciting the warps. In both cases the warp excitation mechanism must be 
ubiquitous. The second point is the warp amplitude. For instance the Milky 
Way has a prominent warp and its amplitude is 2 kpc at Galactocentric dis- 
tances i? ~ 14 kpc(Henderson et al. 1982), and 4 kpc at i? ~ 20 kpc (Diplas 
& Savage 1991). The torque that can bend the disk in such a large amplitudes 
is so strong, which make it difficult to find the source of the torque (Hunter & 
Toomre 1969). The third point is alignment of the warp nodes. Most of spiral 
galaxies including our own have straight line of nodes within the Holmberg 
radius, and outward the lines of nodes make leading spirals (Briggs 1990). 
This fact presents a contrast to the spiral arms in the grand-design spiral. 

Simple analysis (Binney 1992) shows that the simplest interpretation of the 
warp as a vertical oscillation of matters without considering selfgravity, is not 
working, because the line of nodes wind up in a time scale of 1 or 2 Gyr. There- 
fore self-consistent treatment is necessary. There are several hopeful theories, 
all of which have, however, some shortcomings. 

The first step to understand long-lived warps is the normal mode analysis in 
steady state galaxies. Hunter & Toomre (1969) showed that there are many 
bending normal modes in their thin-disk models, but unless density of the 
disks are cut-off sharply at their edges the frequency spectrum is at least 
partly continuous. In such disks any realistic perturbation would excite prop- 
agating wave packets, thus the warping wave is dispersed. A long-lived warp 
can correspond only to a discrete mode. If the disk is truncated at some radius, 
there is only one discrete normal mode, which is called a modified tilt mode 
(Sparke & Casertano 1988). Sparke & Casertano (1988) showed that if disk 
is embedded in a flattened halo there is at most one discrete mode. Shapes 
of the mode show good agreement with observed warps, and even though the 
disk is not in the right shape of the modified tilt mode initially, it dynamically 
evolves to fit the mode (Hofner & Sparke 1994). 
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The normal mode theory works very well as long as the halo is treated as a rigid 
body. However, once the self-consistent response of halo is included several 
problems emerge. For example, dynamical friction acting on the disk from 
the halo dissipates the warping motion (Nelson & Tremaine 1995, Dubinski 
& Kuijken 1995), and also the hne of nodes often winds up tightly (Binney 
et al. 1998, Ideta et al. 2000). To find the normal modes in fully self-consistent 
composite systems of a disk and a halo is extremely intractable. 

One possible mechanism to keep the misalignment between a disk and a halo 
symmetric plane is continuous addition of matter with misaligned angular 
momentum to the halo. Inclined cosmic infall (Ostriker & Binney 1989, Jiang 
& Binney 1999) might support the mechanism. Jiang & Binney (1999) have 
shown that when infall reorientates the outer part of a galactic halo by several 
degrees per Gyr, a realistic warp is excited. Though the result seems appeal- 
ing, their model of the infall is still too idealistic, and numerical studies are 
relatively short (~ 1 Gyr) to explain long-lived warps. Refined studies are 
necessary. 

Another mechanism, which we are examining in this paper, is interaction with 
satellites. This mechanism is intuitively plausible because the most frequent 
m — 1 warp is easily associated with a tidal field from a satellite. Further- 
more increasing observational evidence suggests that most of the galaxies have 
satellite dwarfs. In fact, recently a faint dwarf companion has been discovered 
around NGC 5907, which was considered as a typical example of a warping 
galaxy without interacting companions (Shang et al. 1998). Also statistics 
show positive correlation between warping and existence of interacting com- 
panions (Reshetnikov & Gombes 1998) 

The only but severe problem of the interaction scenario is that the tidal forces 
from satelhtes are usually too small. For example. Hunter & Toomre (1969) 
showed that the direct tidal force from the Large Magellanic Gloud (LMG) of 
10^°Mq at 50 kpc can excite a warp with an amplitude at largest 100 pc at 
16 kpc, while the real Galactic warp is about 3 kpc. 

Against the problem, a remedy was proposed by Weinberg (1998) (hereafter 
W1998). Main point of his idea is to incorporate the halo response to the 
orbiting motion of a satellite into the source of the tidal forces. As the satellite 
moves, rotating density waves which have resonant frequencies to the satellite 
rotation are excited. In a nearly isothermal halo the most dominant wake 
appears at the position of 2:1 resonance, which is at about half distance to the 
satellite orbit. And the mass of the wake is about the same as the sateUite, 
thus its tidal force is much larger than those from the satellite. W1998 showed 
with linear analysis for infinitely thin disks, that the tidal forces from the 
Large Magellanic Gloud (LMG) and its wake are large enough to excite the 
observable warp. 
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In order to prove this idea, fully self-consistent numerical simulations of three- 
dimensional disks are necessary. This is, in fact, very difficult and expensive, 
because we have to treat the small vertical motions of a very thin disk em- 
bedded in a very massive halo, which is extended well beyond the LMC orbit. 
W1998 estimated that 1 miUion particles are necessary to reduce numerical 
noise in order to distinguish meaningful dynamics of the disk. 

This requirement is here accomplished by using a hybrid algorithm, TREE- 
SCF (Vine & Sigurdsson 1998), which particularly suits our problem of warp- 
ing disk dynamics. Furthermore, we make the models as close to the real Milky 
Way - LMC system as possible, so that we can make precise evaluation of the 
warp amplitudes. 

In section 2 we first describe by a brief summary the method to construct the 
Galaxy models using Kuijken & Dubinski (1995), then give some parameters 
of the models which we use in this study. The precise input parameters in 
the Kuijken & Dubinski (1995) codes are listed in the appendix. Section 3 
explains algorithms and parameters of the TREE-SCF code. The results of the 
numerical simulations are given in Section 4. Section 5 is devoted to conclusion 
and discussions. 



2 Galaxy Models 

The aim of this study is to examine the influence of the LMC on excitation 
of the Galactic warp. Therefore we must first construct equilibrium Galaxy 
models that have no warp, and then put a satellite in the LMC orbit. 

To construct equilibrium models for a multi-component system like the Milky 
Way is actually a hard task. Many galactic models are constructed with the 
assumptions that the halo and bulge are static backgrounds, or assuming 
that velocity distribution is Gaussian and solving only the Jeans equations. 
The deviation of these approximations from a true equilibrium causes initial 
transient behavior, in other words relaxation into the true equilibrium. That 
behavior changes the initial distribution especially in the disk. In the disk 
warping problem it is necessary to avoid any disturbance on the disk other 
than that from the satellite. 

A suitable initial model is given by Kuijken & Dubinski (1995) (hereafter 
KD1995), which is nearly an exact solution of the coUisionless Boltzmann and 
Poisson's equations. In their method, the configurations of a halo, bulge, and 
disk are given by means of distribution functions which depend only on inte- 
grals of motion. By the Jeans theorem such distribution functions are already 
steady state solutions of the coUisionless Boltzmann equation. Densities can 
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be calculated from the distribution functions, which depend on the positions 
only through the potential function. Therefore the main task is to determine 
the functional form of the potential by solving Poisson's equation numerically. 



The halo distribution function takes the form 



fhalo{E, L1) 



= < 




(1) 



where is the specific angular momentum about the axis of symmetry, and 
E is the relative energy, which is defined so that = at the edge of the 
distribution (where p — 0)(Binney & Tremaine 1987). This distribution func- 
tion has five free parameters: the potential at the center which appears 
implicitly through the definition of the relative energy, the velocity scale (Tq, 
and three factors A, 5, and C . This distribution is based on Evans's model, 
which has axisymmetric logarithmic potential, but is truncated at a finite ra- 
dius by lowering the distribution function in the same manner as the lowered 
isothermal models. The factors A and B control the system fiattening (g) and 
the core radius (i?c), respectively, and all the three factors are scaled by the 
density scale (pi). 

The bulge distribution function is the same as a King model (Binney & 
Tremaine 1987) and has the form 



This distribution has 3 parameters, is the cutoff potential, pb the central 
bulge density, and iTb the bulge velocity dispersion. 

For these two distributions, the densities are given by analytic functions of it! 
and (Kuijken & Dubinski 1995). 

The disk distribution must depend on 3 integrals of motion, in order to keep 
the triaxial velocity ellipsoid in the disk (Dehnen & Binney 1998). The third 
integral is, however, not known as any analytical function, so that distribution 
function cannot be written as an analytical function. For the sake of conve- 
nience, an approximated distribution function is employed which depends on 
the energy in the vertical oscillations, E^ = ^'(-R, z) — ^'(i?, 0) -|- ^v^, and that 
in planer motion, Ep = E — E^, as well as Lz- E^ is quite well conserved for 



Pb(27r<7^)-3/2exp[(*o - *c)/<t^] {exp[-(£; - *e)/<T^] - 1} 



/bulge (-S) — < 



if £; < 



(2) 







otherwise. 
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stars in nearly circular orbits. The form is 



fdisk{Ep, L^, Ez 



(27r3)V2«;(i?^)aV(i?c)a,(i?c 



cxp 



Ep — Ec{Rc 



E., 



(3) 



where R^ and E^. are the radius and energy of a circular orbit with angular 
momentum L^., and Q and k are the circular and epicyclic frequencies at 
radius R^.. The 'tilde' functions pd, (^r, and are free functions. Therefore 
the density and the radial velocity dispersion are conveniently selected as 



Pdisk(-R, z) 



e-^/^<^erfc i exp 



\^/25R, 



■out , 



-4.6187 



^,{R, 3zd) 



,(4) 



and 

al = aUQ)eM-R/Rd). (5) 



Pd and cTz are iteratively adjusted so that the densities on the mid-plane and 
at height z = Sz^ will agree with those given by eq. (4). The distribution of 
the disk has 6 free parameters: the disk mass, Md, the radial scale length R^, 
the vertical scale height z^, the disk truncation radius Rout, the truncation 
width 5 Rout, and the central velocity dispersion of the disk (Tij(O). 

Kuijken & Dubinski (1995) also give a numerical program to calculate the 
distribution functions and densities. Their algorithm consists of the following 
4 steps, (i) A test potential function ^^^\R,z) is given, (ii) The potential is 
substituted into the density functions which are calculated from the distribu- 
tion functions, (iii) Poisson's equation is solved by using Prendergast & Tomer 
(1970) 's method to determine a new potential functions, (iv) The procedure 
2 and 3 are repeated until the potential functions have converged. 

By using Kuijken & Dubinski (1995)'s method, we construct the Galaxy mod- 
els which satisfy observational properties. We apply the following constraints; 

• solar radius : i?o = 8 kpc 

• circular velocity of the disk at the solar radius : Vc = 220 km/s 

• total surface density within 1.1 kpc of the disk plane : 

Ei.i(i?o) = 71 ± 6M0PC-2 (Kuijken & Gilmore 1991) 

• contribution of the disk material to Ei.i : 

Ed(i?o) = 48 ± 9M0PC-2 (Kuijken & Gilmore 1991) 

• total Galaxy mass within 50 kpc : 

Mtot(< 50kpc) = 5Atli X IO^^Mq (Wilkinson & Evans 1999) 

• total Galaxy mass within 170 kpc : 

Mtot(< 170kpc) = 1.91?:? X IO^^Mq (Wilkinson & Evans 1999) 
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• the halo is spherical : q — 1.0, i.e., A — Q (Ibata et al. 2001) 

Moreover, we assume that 

• disk mass is 5 x IO-^^Mq, 

• disk scale length is 3.2 - 3.5 kpc, 

• disk scale height is 200 - 250 pc, 

• disk edge is 25 - 28 kpc, 

• bulge mass is 15% of that of the disk, 

• and bulge size is about 2 kpc. 
Tabic 1 



Sizes of the components in the Galaxy models, which are used in this paper 





Model L 




Model S 






mass (Mq) 


radius 


mass [Mq) 


radius 


Disk 


5 X 10^0 


3.2kpc 


5 X 10^° 


3.5kpc 


Scale height 




224pc 




245pc 


Disk edge 




25.6kpc 




28kpc 


Bulge 


0.76 X 10^° 


2.21kpc 


0.75 X 10^° 


2.38kpc 


Halo 


5.11 X 10^2 


1539kpc 


8.59 X 10^1 


262kpc 



The halo mass and extent is very important, because in Weinberg's scenario, 
the halo plays an important role as a mediator between a satellite and a disk. 
Observational estimate of the halo mass has, however, a huge error, so that 
we examine two extreme cases: one with the largest halo (Model L), and the 
other with the smallest possible halo (Model S). Construction of these models 
is made by KD1995's code, which is available in their web site, and the input 
parameters are listed up in Appendix A. Some characteristic parameters for 
the two models are given in Tabic 1. The disk scale height and the edge radius 
arc similar in both models. The disk surface densities for the model L and S are 
Sd(-Ro) — 45.3 and 45.5 Mgpc"^, while the total surface densities including 
halo contribution Ei.i(i?o) = 65.0 and 69.8 Mgpc"^, respectively. The halo 
masses for the model L and S arc Aftot(< 50kpc) = 5.58 and 4.93 xlO^^Af©, 
and Mtot(< 170kpc) = 2.1 and 0.9 xIO^'^^Mq, respectively. The determination 
of the parameters was done by manual searching of the input parameters. 
Those are not unique solutions, because the parameter space of the KD1995's 
model is larger than the number of the constraints. The model L and S are 
the first best fitted distributions. The bulge and halo radii listed in the table 
are the tidal radii, where their densities fall to zero. For model L, such a large 
tidal radius is necessary to ensure that the halo profile is nearly isothermal 
beyond the LMC orbit (~ 100 kpc). 

Figures 1 and 2 show the circular velocity profiles of the Model L and S, re- 
spectively. Within the disk extent, the circular velocity profiles are almost the 
same in the both models, but in outer regions where the LMC is orbiting, the 
model L still maintains flat rotation, whereas the model S is nearly Keplerian. 
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Fig. 1. Circular velocity profiles for the model L. Panel (a) shows a small scale profile, 
and (b) shows a larger scale profile. The cross sign (+) shows the observational 
circular velocity at the solar radius. 
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Fig. 2. Circular velocity profiles for the model S. Panel (a) shows a small scale profile, 
and (b) shows a larger scale profile. The cross sign (+) shows the observational 
circular velocity at the solar radius. 

The velocity dispersion profiles of the disk for the model L and S are shown in 
Fig. 3 (a) and (b), respectively. Both profiles are almost the same. The model 
velocity dispersions at the solar radius, (ctr, erg, a^) ~ (41,29,13) km/s, are 
closer to those of the thin disk, and Toomre's Q value has the minimum (1.7 
for the model L, and 1.86 for the model S) around the solar radius. 

Dynamical and numerical stability of the models are shown in Sec. 4.1. 

We treat the LMC rather simply. It is modeled as an extended single particle 
with a spherical Hernquist profile. 



p(r) 



Mlmc^o 
27rr(r + tq)^ ' 



(6) 



where the scale length is chosen to tq = 3.2 kpc, then the half mass radius is 
7.7 kpc. The mass of the LMC is set to Mlmc = 10.6 x lO^M©, which is twice 
of the present LMC mass estimated by Alves & Nelson (2000). This large mass 
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Fig. 3. Velocity dispersion profiles of the disk. Panel (a) shows the model L and (b) 
the model S. Three curves correspond to the radial (R), azimuthal (0), and vertical 
(z) components. 

includes contribution from the SMC, which is about 10% of the LMC, and the 
mass which is possibly stripped by tidal field from the Milky Way, during a 
period of 6 Gyr. 



We need to trace the evolution of the Milky Way - LMC system, at least 
6 Gyr to examine the warp dynamics, but it is very difficult to predict the 
orbital parameters of the LMC 6 Gyr ago. Murai & Fujimoto (1980) have 
made a parameter survey with simple assumptions and found that the LMC 
has been revolving around the Milky Way in an eccentric orbit accompanied 
by the SMC. The present position of the LMC is nearly at a pericenter. With 
its kinematical data (Kroupa & Bastian 1997), the orbital plane of the LMC is 
determined. As a realistic assumption of the LMC orbit, we put the satellite in 
the present orbital plane at a distance of 112 kpc for the model L and 97 kpc 
for the model S. Only tangential velocity is given so that the initial position 
will be at apocenter and the eccentricity is 0.34 and 0.24 for the model L and 
S, respectively. As shown in Fig. 11, the satellite sinks to below 50 kpc at a 
pericenter after 6 Gyr. Since our halo is nearly spherical, the satellite remains 
in the initial orbital plane. 



3 Numerical Methods 



3.1 A hybrid N -body code 



Numerical simulations for our problem require very large computational power 
because of its huge dynamical ranges. One important requirement is to keep 
the disk thin enough to distinguish subtle disk vertical dynamics. Artificial 
2-body relaxation is the main source of harmful, i.e. undesired and unphysical 
heating. In order to reduce this effect a large number of particles is necessary. 
In particular, the relaxation among the disk particles is very large, because 
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the velocity dispersions are small. Simple application of Chandrasekhar's re- 
laxation formula (Binney & Tremaine 1987) reads 

_ 4.2 X lO^Vears ( a \ ' ( nTx 



where A is the Coulomb logarithm, and a, m, and n are velocity dispersion, 
mass, and number density of stars. The critical particle number in numerical 
simulations is 10^, with which the relaxation time is ~ 10^*^ years. 

For the halo, internal relaxation is no problem, because it is much hotter 
than the disk. However, a problem is that the halo is 100 times more massive 
than the disk in the model L. Even small noises of the halo could cause big 
influence on the disk or the satellite. To reduce the noise, a large particle 
number is necessary. 

Another difficulty in our simulations is the huge difference in spatial scales and 
also in dynamics between the disk and the halo. The most important feature 
in the disk is delicate vertical motion, while the most interesting features in 
the halo are large scale wakes, which are excited by the satellite motion. It is 
of course possible to treat the system with a single scheme A'^-body code, like 
a tree-code, but another possibility is to deal with the disk and the halo by 
separate algorithms. 

A hybrid A"-body code, which seems suitable for our system, is the SCF-TREE 
code (Vine & Sigurdsson 1998). The disk and bulge are represented by tree 
particles (e.g., Barnes & Hut (1986)), while the halo is treated by a method 
expanding the potential using bi-orthogonal polynomial series, which is widely 
known as SCF (Self-Consistent Field) method (Clutton-Brock 1973, Hernquist 
& Ostriker 1992, Saha 1993). 

A tree code is a rather popular method for general A"-body problems, because 
of its flexibility of system configuration and relatively high force resolution, 
which is controlled by softening lengths. The whole computational space is 
divided into nested cubic cells to form a hierarchical tree structure. Interaction 
between close particles is calculated by direct summation, but contribution 
from distant particles is calculated only as a group of particles contained in a 
large cell. This grouping is made with a simple criterion. 



where s and d are the size and the distance of the grouped cell. is a control 
parameter, which is known as the tolerance parameter. With this grouping the 
number of calculations reduces to 0{N\ogN) (e.g., Barnes & Hut (1986)). 
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A SCF code is, on the other hand, an expansion method for calculation of 
gravitational forces. In the code density and potential are represented by a 
bi-orthogonal set of basis functions. 

p{r) ^J2AnimPnim(r), and ^(r) ^ ^ Anlm^nlm(r), (9) 
nlm nlm 

where each pair pnim and ^nim satisfies Poisson's equation 

V'^^nlm = ^T^Gpnlm, (10) 

and bi-orthogonal relations 

y Pnlm{r)[^n'l'm'{r)]*dr ^ InlSnn'^ll'^mm', (H) 

where Ini are normalization factors. 

The basic procedure of the SCF method is that first N particles are distributed 
to sample the density, then the coefficients Anim are calculated by expanding 
the density distribution. The potential is determined by the latter equation of 
(9), thus forces can be calculated at any position. The particle positions are 
thus advanced to the next time step. 

The major advantage of the SCF method is that the speed of the computa- 
tion scales as 0{N), which is a smaller burden in increasing particle number 
than those of the direct summation method {0{N^)) or usual tree methods 
(C»(7VlogAr)). 

A fundamental property of the SCF, which is particular among other expan- 
sion methods, is that each particle contributes not as a local gravitational 
source (that is so in other expansion methods), but on global modes of den- 
sity distribution, so that there is no particle-particle interaction in the method. 
On the other words it is more suitable for study of the mode interactions. In 
our problem the most important dynamics in the halo is dynamical friction 
on the sateUite. The classical picture of dynamical friction is due to scat- 
tering of field particles as a 2-body problem (Chandrasekhar 1943, Binney 
& Tremaine 1987). However, modern understanding of dynamical friction is 
rather due to interaction between the massive object and wave modes in the 
field (Tremaine & Weinberg 1984). Therefore it seems reasonable to treat the 
halo by the SCF method. 

In practical computation, the expanded series in eq. (9) should be truncated 
in a small number of terms. Since computational time is proportional to the 
number of terms, it is the most important for the SCF method to find a good 
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basis functions which can fit the system distribution well with the lowest order 
term. There are several basis sets introduced in the literature (Clutton-Brock 
1973, Hernquist & Ostriker 1992, Saha 1993, Syer 1995, Earn 1996, Robijn & 
Earn 1996, Zhao 1996), but it is still terribly difficult to deal with a complex 
system like a disk galaxy. 

The SCF-TREE code has strategic advantages more than compromise in our 
problem. By using a tree method the disk and bulge enjoy high resolution and 
flexibihty to possibly large changes in their shape. We do not need so high 
resolution in the halo, but large number of particles are necessary to reduce 
noise, and the wave modes which are excited by the satellite motion should be 
solved correctly. Both requirements are perfectly accomplished by using the 
SCF method. 

3.2 Parameter Specification 

Our numerical code is fundamentally the same as Vine & Sigurdsson (1998). 
The tree part is based on Hernquist (1987) 's algorithm, and for the SCF part 
Hernquist & Ostriker (1992, hereafter HO)'s basis set is incorporated. There 
are 7 steps in the calculations at each time-step: 

(1) Calculate the coefficients Anim of the terms in the SCF expansion. 

(2) Build the tree from all the tree particles. 

(3) Calculate accelerations of TREE and SCF particles by the SCF expan- 
sion. 

(4) Form the tree interaction lists among TREE particles, and calculate ac- 
celerations of TREE particles by the tree system. 

(5) Form the tree interaction lists between SCF particles and TREE particles, 
and calculate accelerations of SCF particles by the tree system. 

(6) Calculate gravitational force from the satellite on all the TREE and SCF 
particles. Acceleration of the satellite is simultaneously calculated. 

(7) Update positions and velocities of all particles and the satellite. 

We use a common time step for all particles, and the time integration scheme 
is the Leap-Frog. The time steps are fixed to 0.75 Myr and 0.875 Myr for the 
model L and S, respectively. These time steps are one fourth of the central 
free fall time, Tq = \\JS7r/ (2Gpo) ~ 3Myr, which is the shortest time scale in 
the system. The time resolution at the center is marginal, but in most part 
of the disk the time step is an order of magnitude smaller than the vertical 
crossing time of the disk. 

In the TREE part, 2^^(= 512k) particles are used in total, and 448k and 64k 
particles are assigned to the disk and bulge, respectively. They have the same 
mass among each component, and mass ratio between the disk and bulge 
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particles is mbuige/^^'disk = 1-05. All particles have a single fixed softening 
length ^ 50 pc, which is 1/8 of the disk initial thickness. The tolerance 
parameter is set to 6' = 0.7, which is applied to TREE-TREE and also TREE- 
SCF interaction lists. In force calculation, up to quadrupole moments are 
included. 

For the SCF calculations, choice of basis set is the most important. We employ 
the HO basis set. Justification of this choice is shown below. 

The angular dependence of Pnim ^-nd ^nim is expanded in the spherical har- 
monics. 



Pnim{r;r) 



where 



A7rpniir)Yim{e, 0) 



(12) 
(13) 



2l + l{l- 


m 


)! 


47r (1 + 


m\)\ 



;P^I(cos^)e''"'^ 



X < 



-1)"^ (m > 0) 
(m < 0) 



, (14) 



and the radial functions are expanded by using the Gegenbauer (Ultraspheri- 
cal) polynomials. 



Pni{r) 
where 



4n{n + 21 + 2) + {21 + 1) {21 + 3) f ' 



A.'Krl 



(1 + ^2y+5/2 



(1 +py+i/2 



~ ' ^ ~ i -T' 
To r"^ + 1 



(15) 
(16) 

(17) 



and the convention G = 1 is adopted. The pair of zeroth order terms in this 
expansion has the form 

11 1 
Pooo = — TT— -37^, '^•ooo = (18) 
2T:r[l + ry 1 + r 



which is known as the Hernquist profile (Hernquist 1990). This profile is intro- 
duced to fit the light curves of elliptical galaxies. It has a density cusp at the 
center, while there is no cusp in our halo model. This discrepancy causes bad 
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fitting of tlie basis set in the central region. Except for that, the expansion 
with relatively small number of terms reproduces the original distribution very 
well. Figure 4 and 5 show the results of the expansion with the number of the 
radial terms 8 for the model L, and 9 for the model S, respectively. The scale 
lengths of the basis set are ro = 115 kpc for the model L, and ro = 21 kpc for 
the model S, respectively. In both figures the red solid lines show the original 
volume density profiles of the models averaged over spherical shells, and the 
blue dashed lines show those reproduced by truncated expansion of the HO ba- 
sis set. Even with this small number of terms, the reproduced densities fit the 
originals fairly well from r — 0.3 kpc to nearly the end of the halo extent. The 
largest discrepancy appears at the center, where the halo component gives the 
smallest contribution. The density of the expanded series becomes comparable 
to that of the bulge only within 10 pc from the center. We believe that this 
difference has no effect on the warping motion at it! > 10 kpc. In fact we found 
no noticeable structural change at the center in our test simulations owing to 
the difference of the SCF densities from the original equilibrium distribution. 

Other than the HO basis set, there are several different basis sets proposed 
for the purpose of galactic dynamics, including ones with a finite central core 
in density profiles. The most well-known example of the 'cored' profile is that 
based on the Plummer distribution. In this profile, however, density falls very 
quickly with increasing radius (oc r~^), therefore the Plummer basis set gives 
much worse fitting to our halo profile. In addition, Zhao (1996) gives a more 
general class of the basis sets, which include HO and the Plummer basis sets, 
but the only profile that has a central core is the Plummer profile. Therefore 
the HO basis set seems the best within our knowledge. 

Even though the unperturbed halo distribution is approximated well by small 
number of the radial functions, the number of terms that is indeed necessary in 
our simulations is determined so that the density perturbations induced by the 
LMC should be also well approximated. We have made several test calculations 
with increasing the radial and angular functions with the presence of the LMC, 
which is reported in the next section. We judge the suitable number of terms 
by examining the satellite sinking, and select the number of the radial and 
angular functions, n — 8 and I — 12, respectively. Including more terms does 
not affect the satellite sinking history. 

The center of the expansion of the SCF basis set is adjusted to the center of 
mass of the most tightly bound particles. By using an energy criterion, about 
one third of the total halo particles are used to determine the center of the 
expansion. 

The number of particles that is used to sample the halo density is 2^^. There- 
fore we use in total 2^° particles in the Galaxy in addition to one extended 
particle for the satellite. All simulations are made on Pentium III - Linux 
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Fig. 4. Spherically averaged density profile for the halo of the model L, and its 
SCF expansion with the scale length of the basis set tq = 115 kpc. Expansion is 
truncated after n = 8 term. 
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Fig. 5. Spherically averaged density profile for the halo of the model S, and its SCF 
expansion with the scale length of the basis set ro = 21 kpc. Expansion is truncated 
after n = 9 term. 



workstations (800MHz). Typical calculation time is 260 seconds per time step 
and in total about 500 hours to complete the evolution up to 6 Gyr. More 
than 60% of the CPU time is used in the TREE part. 



4 Simulations 
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4-1 Stability of the Galaxy models 



First we examine the stability of our Galaxy models. Without a sateUite, 
these models must be stable and exhibit no warping motion. Moreover we 
have to check the validity of our numerical code. Thus the Galaxy models 
are constructed as described in Sec. 2 with in total 2"^^ particles, and their 
evolution is calculated with the SCF-TREE. 

In order to characterize the disk evolution, we make a modal analysis both of 
the disk surface density and of the vertical distribution of the particles. For 
the disk surface density, 0), angular dependence is decomposed into the 
Fourier series, 

oo 

E(i?, 0) = Eo(i?) + Yl ^rn{R) COs(m|0 - lPm{R)\)- (19) 
m=l 

Numerically, the disk particles are binned into annuli, R e [R^^~^\ R^''^], where 
< < < . . . < R^''^ <...< R^''\ (20) 

then averaged coefficients in the annuli are calculated as follows; 

i)<Ri<i?('=) 

Y miCos{m4>i), (21) 

i)<iJi<iJ('=) 

J2 mism{m(j)i), 

i)<Ri<R('=) 

where AS''^''') = 7r{(7?('''))^ — (i?(^~^))^} is the area of the A;th annulus, and 
rrii, Ri and 4>i are the mass, radius in the disk plane and azimuth of the i-th 
particle. 

Figures 6 and 7 show the surface density modes for the model L and S, re- 
spectively. The panel (a) shows the coefficient Sm(-R) at the initial time. Each 
annulus contains 16,384 particles. The m = profiles show the exponential 
distribution of the disk surface density. The higher modes exist only as noises 
at the beginning. After about 1 Gyr (panel (b) and (c)) the profiles show 
that a bar formed at the center. The bar in the model S is strong and large, 
while that in the model L is 2 or 3 times smaller than in the model S. The 
formation of a bar is expected for models with smaller bulges (Ostriker & 
Peebles 1973, Sellwood & Evans 2001). We would take the realistic value for 
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Fig. 6. Modal analysis of the disk surface density profile for the model L. The panel 

(a) shows the coefficients of the modes m = 0, 2, 3, 4 at the initial time. The panel 

(b) is the same as (a) but after 1.2 Gyr. The panel (c) shows positions of the peaks 
of the modes m > 1 at the same epoch as (b). The symbols denote the same modes 
as those in the panel (a). Each annulus contains 16,384 particles. 

the bulge mass rather than take a larger bulge to prevent the bar instability. 
In fact, the Milky Way does have a bar. Therefore our models correctly reflect 
the same dynamics as the real Milky Way. 

Disk vertical displacement within a annulus R^''~^^ < R < R^''^ can be also 
expanded in a Fourier series; 

4alp(0) = E €^ cos(m|0 - ^L'^l) = E K'^ cos(m0) + 6^'^ sin(m0)] .(22) 



A way to determine the coefficients is least squares fitting, in which the square 
deviation. 
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Fig. 7. Modal analysis of the disk surface density profile for th model S. The nota- 
tions are the same as Fig. 6, but the panel (b) and (c) are at T = 1.4 Gyr. 

becomes the minimum. In this expression, the average, (...), is normalized by 
the local surface density in order to avoid virtual warping modes owing to the 
bar or spiral arms. For the sake of simplicity, we will omit the suffix ^''^ from 
now on. The conditions which minimize 5^ are 



85^ 
daQ 

dam 
dbm 



2j2[(^m'{cos{m'(j)i)) + bm'{sm{m'(j)i))] - 2{zi) = 0, 



(24) 



■ 2 ^ [am' (cos(m0j) cos(m'0j)) + bm' (cos(m0j) sin(m'0i))] — 2{zi cos(m0j)) 

m' 

0, (25) 

■ 2 [am' (sin(m0j) cos(m'0i)) + bm' (sin(m0j) sin(m'0i))] - 2{zi sin(m0i)) 



0. 



(26) 



In the analytical limit where A^-body sampling is fairly well, the following 
relations stand 
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(cos(m'0i) cos(m0i)) 

(sin(m'0j) sm{m4>i)) 
(sin (m'0i) cos(m0i)) = 0. 



for m 

1 for m 
1 for m 

for m' 
for m' 



1 

k 2 



7^ m 

= m 7^ 

= m = 0, 

7^ m 
— m, 



(27) 

(28) 
(29) 



If this is the case, the coefficients are determined by the following equations 



ao = {zi), am = 2{ziCos{m(f)i)), b„i = 2 {zi sm{m4>i)) 
(5^ = ^ (^Zi - cos(m0) + bm sin(m0))) 



(30) 



(31) 



Here S is considered as the disk thickness, which is corrected by subtracting 
the effect of warping. 

In practical calculations, however, finite A'"-body sampling causes errors in eqs. 
(27) - (29), the least mean square conditions eqs. (24) - (26) become thus a 
matrix equation. In order to estimate the errors owing to the coarse sampling, 
we also find another form of approximate solutions, in which the higher order 
coefficients are progressively determined by the lower order terms. 



do = 



(32) 



(2;jcos(0j)) - ao(cos(0j)) (zj sin(0i)) - ao(sin(0j)) 



(Cos2(0j)) 



(sin^(0O> 



m—l 



{ziCOs{m(j)i)) - [arn'{cos{m(j)i)cos{m'(f)i)) + brn'{cos{m(j)i)sm{m'(f)i))] 

m'=0 

(cos^(m0i)) 

m—l 

{ziSm{m(f)i)) - ^ [a^/(sin(m0j)cos(mVj)) +6^/(sin(m0j)sin(m'0j))] 



(34) 



m'=0 



(sin (m0j)) 

Differences between eq. (30) and eqs. (32) - (35) are considered as errors. 
Figure 8 shows the amphtudes of the warping mode, hm{R) = ^Ja^{R) + b'^{R), 
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Fig. 8. Amplitudes of the warping modes of the disk in the model L (panel a) and 
the model S (panel b), at T = 1.2 Gyr and 1.4 Gyr, respectively. The modes m = 0, 
1, 2 are shown by the curves with plus (+), diamond (o), and circle (o) symbols, 
respectively, and the initial disk thickness is shown by a curve with box (□) symbols. 
The error bars denote difference between the solutions by eq. (30) and eqs. (32) - 
(35). 
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for m = 0, 1, 2, after about 1 Gyr evolution without a satellite. The errors 
are shown by the error bars. Large errors occur at i? ~ 3 kpc in the model 
S, where the strong bar appears. In other regions the errors are considerably 
smaller, and then invisible in the figures in those scales. 

These warping modes are excited only by discreteness noises. W1998 suggested 
that the numerical noises could raise a large warp, if the number of particle 
in the system is less than 1,000,000. Comparing with the initial disk heights 
shown with light blue curves, it is obvious that the model galaxy disks are 
kept sufficiently flat over 1 Gyr. 

Two-body relaxation within the disk, and also random noise owing to the 
discrete distribution of halo material cause heating of the disk. This effect 
is estimated by change in velocity dispersions of the disk particles. Fig. 9 
shows the evolution of <7z{R), which is evaluated in the same annuli as in 
the warp modal analysis for the model L and S. In both models, o"^ increases 
significantly in the central part. This is explained by the disturbance from the 
bulge. In particular in the model S, the formation of the strong bar and its 
rotating motion might heat up the disk distribution. Except for the central 
part (< 5 kpc), the increase in (Tz is about 10% (or 1% in terms of energy 
increase), after about 1 Gyr evolution. In the main simulation we will follow 
the evolution six times longer, but expected heating would be still harmless. 

Prom these analysis we are convinced that our Galaxy models are stable 
enough to discriminate the disk dynamics due to an orbiting satellite. 

4-2 Satellite sinking and dynamical friction 

The halo distribution is represented by truncated SCF basis-set. The necessary 
number of terms is determined by the requirement that the response of the 
halo to the sateUite motion should be correctly solved. In this subsection we 
examine the nature of the SCF with respect to the number of terms included 
in the basis-set expansion. 

We put the LMC model with a mass of 1.1 x 1O^°M0, initially in a circular 
orbit with a radius it! = 64 kpc in the symmetry plane of the model L, but the 
bulge and disk component are replaced by a fixed external field. Only the halo 
component is solved by the SCF method. The changes in radius according to 
the satellite sinking are shown in Fig. 10. 

I have made six different runs with different number of terms in the expansion. 
As we have seen in Fig. 4, truncation of the radial expansion terms at n = 8 
or 9 yields a very good fit to the unperturbed halo distribution. On the other 
hand, the dynamical friction acting on the satellite should be more sensitive to 
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Fig. 9. Change in the velocity dispersion in z-direction. The panel (a) and (b) show 
vertical velocity dispersion to the disk plane (a^) in the model L and S, respectively. 

the number of the azimuthal expansion terms. Fig. 10 (a) shows the sateUite's 
radius evolution calculated with different / of the spherical harmonics. All the 
three runs are made with the radial expansion terms up to n = 8, but the 
spherical harmonics are included up to / = 6, 12, and 18. For each /, the 
parameter m runs from —I to I. In the run with I = 6 the satellite sinking rate 
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Fig. 10. Evolution of the satellite's orbits calculated with different SCF parameters. 
Each parameter of n and / means that the terms up to that value of n and / are 
included in the expansion of the basis set. 

is slower than those with larger /. This is because contribution of terms with 
higher / to dynamical friction is missing. Increasing / from 12 to 18 makes the 
sinking rate faster, but the difference is not so large in comparison with the run 
with / = 6. From the figure it is surmised that dynamical friction with / = 12 
is not far from the true value. The number of the spherical harmonics up to 
/ = 12 is the best compromise between accuracy of the dynamical friction and 
economy of the computational time. 

In Fig. 10 (b) another comparison regarding the effect of the radial expansion 
number n is shown. We find that including terms larger than n = 8 does 
not change dynamical friction so much. Furthermore, we can make sure that 
number of the azimuthal expansion terms is the most important for dynamical 
friction because the run with n = 16 but / = 6 exhibits still smaller sinking 
rate. 

4-3 Warp excitation 

Finally we make simulations of the fully live Milky Way models interacting 
with the LMC model. Now the satellite is put in the realistic orbit described 
in section 2. The satellite is sinking in the halo owing to dynamical friction. 
The time variations in the distance of the LMC from the center of the Milky 
Way in the model L and S are shown in Fig. 11. In both models, the LMC 
is at a pericenter closer than 40 kpc near 6 Gyr. In fact, the pericenter at 
40 kpc is smaller than that of the real LMC. This choice of the orbit makes 
the gravitational tide a little bit stronger in our model simulation. We would 
consider that the moment of 6 Gyr corresponds to the present day status. One 
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Fig. 11. Change in the radius of the satelhte in the model L (a) and the model S 
(b). 

might notice that in the model L the LMC sinks regularly, whereas it does not 
in the model S. Since the halo in the model S has the tidal cut off at 262 kpc, 
the satellite is on the edge of the halo distribution. This is because the halo 
in the model S has the tidal cut off around the orbit of the LMC. 



Time variation in the disk warping modes in the model L are shown in Fig. 12. 
Three rows (a), (b), and (c) correspond to the time T = 0, 3.0, and 6.0 Gyr. 
In each row, the left panel shows the amplitudes of the warping mode with 
m < 3, together with the disk thickness. The right panel shows the position 
angle of the peak of the warping modes for each annuli in the disk plane. 
The diamond (blue) and round (pink) symbols denote m = 1 and 2 modes, 
respectively. The disk rotation is clockwise in the figures. The LMC orbital 
plane intersects the disk nearly along the y-axis in the right panel. The first 
two intersections of the LMC with the initial disk plane are at (X, Y) = {6.8, 
-91.8) kpc and (-6.8, 89.3) kpc. At the beginning (panel a) all the warp modes 
are only due to noise, and thus the position angles are randomly distributed. 

The observed warp amplitudes are also shown in the left panels for comparison. 
For the Milky Way, not only the m = 1 but also m = and 2 modes are known 
to exist. We take an approximated formula for the warp amplitude given by 
Binney & Merrifield (1998), with a correction for the solar radius Rq = 8 kpc, 

^, i?/kpc-10.4 . ^ ^ ^ /i?/kpc - 10.4 V ,^ /o^N 

^warp,obs(^,0) = ^ Sm 0-^0.3 I — I (1 - COS 20). (36) 



In Fig. 12 and 13 the observed m = 1 and 2 modes are shown by dashed and 
dotted curves, respectively. 

At 3 Gyr (panel b) disk warping appears in several modes. The most noticeable 
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mode is m = (red solid line with "+" symbol). It has a wavy figure with 
the amplitudes about 250 pc upward and 400 pc downward. All the other 
modes have simple shapes with nearly same amplitudes. The amplitudes grow 
outwards, but the biggest amplitudes are still less than 500 pc. The m — 1 
mode shows a clear bending in the amplitude at around R — 10 kpc, and the 
peaks of the m — 1 mode make a broad ridge, which slightly winds towards 
the same direction as the disk rotation. The m = 2 global bending mode 
is clearer in the position angle diagram which starts at 5 kpc, and winds in 
trailing sense. But outside of 11 kpc the mode becomes dispersed again, and 
shows no coherent feature. 

At 6 Gyr (Panel c), which is considered as the present state, the amplitudes 
of the m = and 1 modes become larger. The largest amplitude of the m = 
is about 1 kpc. A bend in the m — 1 curve appears aX, R — 7.5 kpc, 
growing outwards up to the largest amplitude of 1 kpc. The ridge of the m — 1 
mode winds in trailing sense until 15 kpc, then turns to leading sense. Similar 
winding feature appears in the m = 2 mode, which starts at 10 kpc, winds in 
trailing sense until 10 kpc, then starts winding oppositely. Since the LMC is 
revolving retrograde to the disk rotation, I surmise that warping modes are 
affected by the trailing wakes of the LMC (see Fig. 14 to Fig. 17), and that is 
the reason of the change of the trend in winding of the warping modes. 

The most remarkable result is that the warp amplitudes we have got from 
the simulation are already comparable to the observed ones. For the m = 
mode the amplitude is nearly the same amplitude as the observed one, though 
the sign is opposite. This negative amplitude in the m = might depend 
on the initial position of the LMC, which we took arbitrary in the present 
orbital plane. The m = 1 and m = 2 modes have more than a half of the 
observed amphtudes. The observed trend that the m = 2 mode has about 
a half amplitude of the m — 1 mode is also realized in the simulation. This 
result proves that in cooperation with the halo response the LMC can create 
enough tide to excite observable warps. 

It should be noted that the disk plane is getting inclined homogeneously with 
respect to the initial disk plane. The tilting angles of inner part of the disk 
are about 3 degrees at T = 3 Gyr, and 6 degrees at T = 6 Gyr. Since 
the homogeneous tilting is not observed as a warp, these tilts are already 
subtracted from the m = 1 mode. In fact, tilting itself can be the cause of 
warp (Sparke & Casertano 1988) if a halo is flattened. Our halo models are, 
however, spherical, so that tilting cannot develop the integral-shaped warps. 

Fig. 13 shows the warping modes in the model S. The notations are the same 
as those in Fig. 12. At T = 3 Gyr (panel b), a coherent m = 1 bending mode 
(dashed curve with filled diamond symbol) is developed. The bend starts at 
8 kpc, and the ridge of the m — 1 warp makes a straight line. This feature 
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Fig. 12. Warp modes in the model L. The rows (a), (b), (c) correspond to the 
time T = 0, 3.0, and 6.0 Gyr, respectively. In each row, the left panel shows the 
amplitudes of the warping mode, and the right panel shows the position angle of 
the peak of the warping modes in the disk plane. The disk rotation is clockwise in 
the figures. The LMC orbital plane intersects the disk nearly along the y-axis in 
the right panel. The amplitude for uniform tilting of the disk is already subtracted 
from m = 1 amplitudes in these figures. 



resembles the observed warp very well, but its amplitude is at largest 400 pc. 
This is significantly smaller than the observed warp. The other modes are much 
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Fig. 13. Warp modes in the model S. Notations 
Fig. 12. The error bars seen in the left panel (c) are 
correction on the estimate of warp modes, which is 
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in this figure are the same as 
the errors in the surface density 
explained in section 4.1. 
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smaller than the m = 1 mode. At T = 6 Gyr (panel c), all the warping modes 
still have small amplitudes, thus the development of warps seems already 
saturated. In the m = 1 warping mode there is a step in the amplitudes at 
15 kpc. Within this radius the peak position angles align in a straight line, but 
outside of the radius, the peak position angles wind and folded. This might 
be also due to the dominant effect of the wakes in the halo. 
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One thing that is clear by comparing the results of the two models is the 
importance of the halo distribution for excitation of warps. Both models have 
the same halo density at the solar radius, but the halo in the model S has a 
smaller cut off radius, hence the halo in the model S has smaller densities in 
the region where the LMC is orbiting. The wakes excited in the model S are 
thus smaller than those in the model L. 
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Fig. 14. Density contour of the m = 1 wake in the model L, in the section of the 
orbital plane of the LMC. The position of the LMC is denoted by a big dot, where 
its trajectory is shown by the dashed curve. The levels of the contours are shown 
in logarithmic density scales only for positive density enhancement. 
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Fig. 15. Density contour of the m = 2 wake in the model L. Notation is the same 
as Fig. 14. The enhancements that align to the horizontal axes in the central region 
with i? < 20 kpc are contribution from the disks. 

Figs. 14 to 17 show contour maps of the m = 1 and 2 wakes in the section 
of the halo cut by the orbital plane of the LMC. In the figures the disks are 
nearly edge on with the inclination of 80 degrees, and the section of the disk 
planes are aligned to the horizontal axes. The wakes are defined as spherical 
harmonic components of the halo density. Only positive density regions are 
plotted. The position of the LMC is shown by a large dot, and its trajectory 
to 6 Gyr is superposed with a dashed curve. In this plane the LMC revolves 
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Fig. 16. Density contour of the m = 1 wake in the model S. Notation is the same 
as Fig. 14. 
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Fig. 17. Density contour of the m = 2 wake in the model S. Notation is the same as 
Fig. 14. The enhancements that align to the horizontal axes in the central region 
with i? < 20 kpc are contribution from the disks. 



anti-clockwise. 



One should note that the amplitude of the m = 1 mode has uncertainty de- 
pending on the selection of the origin of expansion. Even for the spherical 
distribution, deviation of the origin from the center of the distribution yield 
virtual m = 1 mode. For example, the deviation of the origin of 625 pc causes 
a m = 1 mode with amplitude of lO~^M0pc~^ at r = 50 kpc, which is compa- 
rable to the one shown in Fig. 14. Our choice of the origin of expansion is the 
center of mass of the most tightly bound particles, which is determined in the 
same way as the SCF expansion. With this choice of the origin, the error in 
the determination of the origin of expansion cannot be larger than the size of 
the core (~ 500 pc), so that the errors in the amplitude of the wake are much 
smaller than the amplitude of the wake. 

The innermost density enhancements with a size of 20 kpc, which are seen in 
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the m = 2 plots (Fig. 15 and 17), are contributions from the disk. These figures 
show that large scale wakes are excited nearly at the radius of the LMC. It is 
clear that the wakes in the model L are larger and stronger than those in the 
model S. In particular for the m — 2 wakes, which would be the direct cause 
of the integral-shaped warps, two crescent-shaped wakes are seen in Fig. 15. 
One is an outer pair with a radius of 70 kpc, which typically appears in the 
lO^^Mo/pc'^ contours, and the other is an inner pair with a radius of 30 kpc, 
which are seen as the contours of lO~^M0/pc^. Those should correspond to 
the wakes which have resonant rotational frequencies 1:1 and 2:1 to that of 
the LMC. This is exactly the same feature as those obtained in Weinberg's 
analysis. In the model S (Fig. 17), however, no clear inner wake can be seen 
other than the disk contribution. This might be because the halo density is 
too small and the disk gravity is still dominating in the region. 



5 Conclusion and Discussions 

Our study is meant to be an extension of Weinberg's linear analysis (Weinberg 
1998), in the way that (i) the simulations are made by a fuUy-selfconsistent 
A^-body code, (ii) the disk is three-dimensional, and (iii) the Galaxy models 
are more realistic. The warp amplitude that we have got in the model L is 
nearly the same as that of the largest warp in W1998. An interesting difference 
between our results and W1998 is that our smaller halo mass model yields 
smaller warp amplitudes, while in W1998 larger halo models yield smaller 
warps. This is, in fact, explainable as follows. In W1998, different halo masses 
produce different rotation curves, so that larger halo mass models reduce the 
efficiency of resonant amplification of warps. In our case the masses of the disk, 
bulge and halo within = 50 kpc are nearly the same between the models, 
hence the rotation curves in smaller scales are also the same. Therefore the 
resonant structures are not very different between the models. On the other 
hand, the mass within R = 170 kpc is twice as massive in the model L as that 
in the model S. In fact the warp amplitude in the model L is roughly double 
of the model S. This linear dependence shows that the mass of a halo in the 
region where a sateUite is orbiting has a direct effect on warp amphtudes. As 
a result we confirmed W1998's prediction that halos play crucial role in warp 
excitation. 

Our high resolution simulations upon the Galaxy models with reasonably re- 
alistic parameters make it possible to compare the numerical results with the 
observations seriously. Our larger halo model yields very large warp, which is 
comparable to the observed warp. The m = harmonic mode has nearly the 
same amplitude as the observed one, while those of m = 1 and 2 are about a 
half. These values are in fact remarkably close to those of the observed warp. 
Now the LMC has revived as a possible candidate for the main cause of the 
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Galactic warp. Though there are still problems, e.g., the Unes of nodes are 
not straight, we have resolved the most serious problem of the interaction sce- 
nario, that the gravitational tide from the LMC would be too small to excite 
the observed warp. 

In this paper, we examined the LMC as a main contributer to the Galac- 
tic warp. On the other hand, the Sgr dwarf was posed as another candidate. 
Though its current mass is much smaller (~ IO^Mq) than the LMC, its po- 
sition is closer (20 kpc), thus the gravitational tide is comparable to that of 
the LMC. Jiang & Binney (2000) proposed a possibility that the Sgr dwarf 
had ~ 10^^ initially. Moreover, since the orbital plane of the Sgr dwarf 
is perpendicular to that of the LMC, the induce warp would have a line of 
nodes perpendicular to that induced by the LMC. Therefore examination of 
different satellite parameters is certainly important to find out the origin of 
the Galactic warp. 

We have also found that our smaller halo model cannot bear the observed 
warp. This result shows that even in the interaction scenario the halos play 
an important role, so that we could use the warp kinematics as a probe of 
the halo structure. Our results might suggest that a massive halo is neces- 
sary to excite the observed Galactic halo if its origin is the LMC. In order 
to restrict halo properties, however, we need to explore wider model param- 
eters. We could change the contribution of the halo material within the disk 
radius smaller to examine the maximal disk models. While Mera et al. (1998) 
argued that a maximal disk model is excluded, a possibility that the Milky 
Way has a maximal disk is shown in several papers(Sackett 1997, Palunas & 
Williams 2000, Binney & Evans 2001). Moreover the strongest warp appeared 
in the maximal disk model in W1998, which motivated us to investigate the 
maximal disk models. Another possible extension of model parameters is the 
halo fiattening. Our selection of the spherical halo (g = 1.0) is based on the 
analysis of Ibata et al. (2001), but other analysis shows that the most probable 
flatness of the halo is c/a = 0.8 (Oiling & Merrifleld 2000). As mentioned in 
Sect. 4.3 the disk in the model L is getting inclined to not negligible extent, so 
that if the halo is flattened the tilting disks would get additional torque from 
the halo, which might make the warp amplitude larger. Since the shape of the 
halo is still under big arguments, for example, studies of galaxy formation by 
means of cosmological simulations conclude that the galactic halos are much 
more flattened such as c/a ~ 0.5 (Warren et al. 1992, Dubinski 1994). Not 
only the shape but also the different velocity distribution in the halo may 
affect the warp amplitude, because the wake structure might be different. We 
will tackle these problems in detail in subsequent papers. 

There is another interesting outcome from our simulations besides the warp. 
Quite large disk heating has taken place especially in the model L. Even though 
we cannot exclude the artifical heating effect in numerical schemes, since such 
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a large heating is not observed in the test simulation without a satellite. We 
surmise that interaction with a satellite also causes disk heating. The same 
heating owing to interaction with a satellite was studied by Velazquez & White 
(1999), who claimed that sateUites on prograde orbits causes disk heating more 
than those on retrograde orbits. Our simulations results larger heating than 
Velazquez & White (1999). This is partly because our initial disks are about 
three times thinner than those in Velazquez & White (1999). Since the disk 
heating is, as mentioned in Sect. 3.1, a very delicate problem for all numerical 
methods, we need much more careful treatment for this problem. This will be 
our next project. 
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A Input parameters for the Galaxy model construction 



For construction of the initial Galaxy models, I have used the software 'Galac- 
tlCS' (Kuijken & Dubinski 1995). This software is available on their web site 
( Jhttp : //www. astro . rug.nl/~-i.>kuijken/galactics .htmlj ) . The software re- 
quires several input parameters. Here is the list of the input parameters. The 
parameters explained in Sec. 2 are listed in the column 1, and the correspond- 
ing variables used in the codes are listed in the column 2. The column 3 and 
4 are the values of the parameters for the model L and S, respectively. All the 
values are in the computational units, where the gravitational constant, the 
mass of the disk, and the scale length of the disk are set to unity. 

Table A.l 

Input parameters given to the software GalactlCS. 







model L 


model S 


halo 








cental potential ^'o 


psiOO 


-5.0 


-4.0 


velocity scale uo 


vO 


0.90 


1.0 


potential flattening q 


q 


1.0 


1.0 


core parameter R'^/R^ *^ 


coreparam 


0.1 


0.1 


characteristic halo radius i?a *^ 


ra 


0.7 


0.8 


disk 








disk mass 


rmdisk 


1.0 


1.0 


scale length R^ 


rdisk 


1.0 


1.0 


disk radius Rout 


outdisk 


7.0 


7.0 


scale height 


zdisk 


0.07 


0.07 


truncation width 5R 


drtrunc 


0.5 


0.5 


central radial vel. dispersion, o"/j(0) 


sigrO 


0.55 


0.55 


scale length of a\ 


disksr 


1.0 


1.0 


bulge 








central density pb 


rhol 


15.0 


15.0 


cutoff potential ^'c 


psiout 


-4.12 


-3.1 


velocity dispersion Ub 


sigbulge 


0.8 


0.8 


radial step size 


dr 


0.01 


0.01 


number of bins 


nr 


51000 


10000 


Maximum azimuthal harmonic 


Imax 


10 


10 



(*1) i^K is the King radius 
(*2) Characteristic radius is related to the density scale pi by the following relation; 
i?a = [3/27rGpi 



;i]VVoe*o/2-g 
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